Loading the data¶
In [1]:
%matplotlib inline
import sys
import re
import numpy
DO_CORRELATIONS, DO_PROJECTIONS, DO_REGRESSIONS, DO_CLUSTERS = True, True, True, True
if DO_PROJECTIONS:
# from sklearn import manifold
from sklearn import decomposition
if DO_REGRESSIONS:
import statsmodels.api as sm
if DO_CLUSTERS:
from sklearn import cluster
from sklearn import metrics
import scipy.cluster
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt
import pdb
from tools_all import *
NBCS = [3,5,7]
NBC = max(NBCS)
cmD = make_colormap([d[1] for d in TOL_COLORS if d[0] <= NBC], extrapolate=False)
cmC = make_colormap([d[1] for d in TOL_COLORS], extrapolate=True)
FSIZE = 3
FIGS_REP = "figs/"
FIGS_EXT = ".png"
# FILE_AGG="data_focus_agg.csv"
# FILE_BIO="data_focus_bio.csv"
FILE_AGG="IUCN_EU_nbspc3+_focus_agg.csv"
FILE_BIO="IUCN_EU_nbspc3+_focus_bio.csv"
FILE_CLU="data_clusters.csv"
groups_bio = [['bioX1:NPP', 'bio15:PSeason', 'bio12:PTotY'],
['bio13:PWetM', 'bio16:PWetQ', 'bio14:PDryM', 'bio17:PDryQ', 'bio19:PColdQ', 'bio18:PWarmQ'],
['bio4:TSeason', 'bio7:TRngY', 'bio2:TMeanRngD', 'bio3:TIso', 'bio1:TMeanY'],
['bio5:TMaxWarmM', 'bio10:TMeanWarmQuarter', 'bio6:TMinColdM', 'bio11:TMeanColdQ', 'bio8:TMeanWetQ', 'bio9:TMeanDryQ']
]
groups_agg = [['MEAN_HYP', 'MEAN_LOP', 'MEAN_OO', 'MEAN_AL', 'MEAN_OL', 'MEAN_SF'],
['MEAN_HOD', 'MEAN_LOPT', 'MEAN_OT', 'MEAN_CM', 'MEAN_ETH']
]
ex_agg = ['MEAN_HOD', 'MEAN_LOPT'] #, 'MEAN_OT', 'MEAN_CM']
# groups_bio = [['bioX1:NPP', 'bio15:PSeason', 'bio12:PTotY']]
# groups_agg = [['MEAN_HYP', 'MEAN_LOP', 'MEAN_OO', 'MEAN_AL', 'MEAN_OL', 'MEAN_SF']]
ABC = load_ABC(FILE_AGG, FILE_BIO, FILE_CLU)
A, head_agg, map_agg = ABC[0]["data"], ABC[0]["head"], ABC[0]["map"]
B, head_bio, map_bio = ABC[1]["data"], ABC[1]["head"], ABC[1]["map"]
C, head_clu, map_clu = ABC[2]["data"], ABC[2]["head"], ABC[2]["map"]
#### Log precipitations
B_log = numpy.vstack([numpy.log10(B[:,i]+1) if re.match("bio[0-9]+:P", hb) else B[:,i] for i, hb in enumerate(head_bio)]).T
# org_map_vnames = dict(map_vnames)
# for i, hb in enumerate(head_bio):
# if re.match("bio[0-9]+:P", hb):
# map_vnames[hb] = "\\log(%s)" % org_map_vnames[hb]
#### filter variables
keep_agg = [vi for vi, v in enumerate(head_agg) if re.match("MEAN_", v)]
keep_aggM = [vi for vi, v in enumerate(head_agg) if re.match("MEAN_", v) and v != "MEAN_HOD" and v != "MEAN_LOPT"]
keep_bio = [vi for vi, v in enumerate(head_bio) if re.match("bio[0-9]+:", v)]
Akw = withen(A[:, keep_agg])
AMkw = withen(A[:, keep_aggM])
Bkw = withen(B[:, keep_bio])
BLkw = withen(B_log[:, keep_bio])
hA = [head_agg[i] for i in keep_agg]
hAM = [head_agg[i] for i in keep_aggM]
hB = [head_bio[i] for i in keep_bio]
# clu_data = BLkw
# clu_data = AMkw #Akw
clu_data = numpy.hstack([BLkw, AMkw])
lng_lat = B[:, [map_bio["longitude"], map_bio["latitude"]]]
In [2]:
def mk_plot(oname):
# plt.savefig(FIGS_REP+oname+FIGS_EXT)
plt.show()
# pass
def mk_out(fname):
# return open(FIGS_REP+fname, "w")
return sys.stdout
def close_out(fo):
# fo.close()
pass
In [3]:
if DO_CORRELATIONS:
Avs = []
Ahs = []
for hs in groups_agg:
Avs.extend([map_agg[ha] for ha in hs])
Ahs.extend(hs)
n, D, vs, hs = ("Dent", A, Avs, Ahs)
corrs = numpy.corrcoef(D[:,vs].T)
nz = len(vs)
fig, axes = plt.subplots(len(vs)-1, nz-1, figsize=(FSIZE*nz, FSIZE*nz),
sharex='col', sharey='row', gridspec_kw={'hspace': 0, 'wspace': 0})
plot_corrs(D, corrs, vs, hs, axes)
mk_plot("correlate%s" % n)
Climate variables¶
In [4]:
if DO_CORRELATIONS:
Bvs = []
Bhs = []
for hs in groups_bio:
Bvs.extend([map_bio[hb] for hb in hs])
Bhs.extend(hs)
n, D, vs, hs = ("Clim", B_log, Bvs, Bhs)
corrs = numpy.corrcoef(D[:,vs].T)
nz = len(vs)
fig, axes = plt.subplots(nz-1, nz-1, figsize=(FSIZE*nz, FSIZE*nz),
sharex='col', sharey='row', gridspec_kw={'hspace': 0, 'wspace': 0})
plot_corrs(D, corrs, vs, hs, axes)
mk_plot("correlate%s" % n)
Projections¶
In [5]:
################################
def project_data(X, model, mname, head, dname, axes, di):
U = numpy.eye(X.shape[1], X.shape[1])
if mname == "PCA":
model.fit(X)
Xproj = model.transform(X)
Uproj = model.transform(U)
else:
Xproj = model.fit_transform(numpy.vstack([X[::50, :], U]))
Uproj = Xproj[-U.shape[0]:, :]
plot_proj(Xproj, Uproj, head, dname, axes, di)
################################
if DO_PROJECTIONS:
NCP=2
datas = [("Dent", Akw, hA),
("Dent-{HOD, LOPT}", AMkw, hAM),
("Clim", BLkw, hB),
("Dent+Clim", numpy.hstack([Akw, BLkw]), hA+hB)
]
models = [("PCA", decomposition.PCA(n_components=NCP))]#,
# ("PCoA", manifold.MDS(n_components=NCP, metric=True)),
# ("MDS", manifold.MDS(n_components=NCP, metric=False))]
for mi, (mname, model) in enumerate(models):
fig, axes = plt.subplots(2, len(datas), figsize=(FSIZE*len(datas), FSIZE*2))
for di, (dname, X, head) in enumerate(datas):
project_data(X, model, mname, head, dname, axes, di)
U = numpy.eye(X.shape[1], X.shape[1])
mk_plot("projection%s" % mname)
In [6]:
################################
def regress_one_groups(gA, gB, A, head_A, map_A, B, head_B, map_B, collect={}, axes=None):
for i, hb in enumerate(gB):
if hb not in collect:
collect[hb] = {}
for j, ha in enumerate(gA):
ai = map_A[ha]
bi = map_B[hb]
##################################
# y = A[:,ai]
# X = B[:,[bi]]
# defv = def_vals.get(ha, 0)
# if defv is not None:
# ### drop rows with agg==0
# X = X[y!=defv, :]
# y = y[y!=defv]
##################################
X = A[:,[ai]]
y = B[:,bi]
Xr = sm.add_constant(X)
# model = sm.OLS(y, Xr)
model = sm.GLS(y, Xr)
results = model.fit()
if axes is not None:
plot_reg_one(X, y, results, ai, head_A, bi, head_B, axes, i, j)
collect[hb][ha] = results
return collect
################################
In [7]:
if DO_REGRESSIONS:
TOPR = 3
collect = {}
gbi, gai = (0, 0)
gB = groups_bio[gbi]
gA = groups_agg[gai]
nx, ny = len(gA), len(gB)
# ny, nx = len(gA), len(gB)
fig, axes = plt.subplots(ny, nx, figsize=(FSIZE*nx, FSIZE*ny),
sharex='col', sharey='row', gridspec_kw={'hspace': 0, 'wspace': 0})
regress_one_groups(gA, gB, A, head_agg, map_agg, B_log, head_bio, map_bio, collect, axes)
mk_plot("regression1VDent%dvClim%d" % (gai, gbi))
In [8]:
if DO_REGRESSIONS:
gbi, gai = (0, 1)
gB = groups_bio[gbi]
gA = groups_agg[gai]
nx, ny = len(gA), len(gB)
# ny, nx = len(gA), len(gB)
fig, axes = plt.subplots(ny, nx, figsize=(FSIZE*nx, FSIZE*ny),
sharex='col', sharey='row', gridspec_kw={'hspace': 0, 'wspace': 0})
regress_one_groups(gA, gB, A, head_agg, map_agg, B_log, head_bio, map_bio, collect, axes)
mk_plot("regression1VDent%dvClim%d" % (gai, gbi))
In [9]:
if DO_REGRESSIONS:
gbi, gai = (1, 0)
gB = groups_bio[gbi]
gA = groups_agg[gai]
nx, ny = len(gA), len(gB)
# ny, nx = len(gA), len(gB)
fig, axes = plt.subplots(ny, nx, figsize=(FSIZE*nx, FSIZE*ny),
sharex='col', sharey='row', gridspec_kw={'hspace': 0, 'wspace': 0})
regress_one_groups(gA, gB, A, head_agg, map_agg, B_log, head_bio, map_bio, collect, axes)
mk_plot("regression1VDent%dvClim%d" % (gai, gbi))
In [10]:
if DO_REGRESSIONS:
gbi, gai = (1, 1)
gB = groups_bio[gbi]
gA = groups_agg[gai]
nx, ny = len(gA), len(gB)
# ny, nx = len(gA), len(gB)
fig, axes = plt.subplots(ny, nx, figsize=(FSIZE*nx, FSIZE*ny),
sharex='col', sharey='row', gridspec_kw={'hspace': 0, 'wspace': 0})
regress_one_groups(gA, gB, A, head_agg, map_agg, B_log, head_bio, map_bio, collect, axes)
mk_plot("regression1VDent%dvClim%d" % (gai, gbi))
In [11]:
if DO_REGRESSIONS:
gbi, gai = (2, 0)
gB = groups_bio[gbi]
gA = groups_agg[gai]
nx, ny = len(gA), len(gB)
# ny, nx = len(gA), len(gB)
fig, axes = plt.subplots(ny, nx, figsize=(FSIZE*nx, FSIZE*ny),
sharex='col', sharey='row', gridspec_kw={'hspace': 0, 'wspace': 0})
regress_one_groups(gA, gB, A, head_agg, map_agg, B_log, head_bio, map_bio, collect, axes)
mk_plot("regression1VDent%dvClim%d" % (gai, gbi))
In [12]:
if DO_REGRESSIONS:
gbi, gai = (2, 1)
gB = groups_bio[gbi]
gA = groups_agg[gai]
nx, ny = len(gA), len(gB)
# ny, nx = len(gA), len(gB)
fig, axes = plt.subplots(ny, nx, figsize=(FSIZE*nx, FSIZE*ny),
sharex='col', sharey='row', gridspec_kw={'hspace': 0, 'wspace': 0})
regress_one_groups(gA, gB, A, head_agg, map_agg, B_log, head_bio, map_bio, collect, axes)
mk_plot("regression1VDent%dvClim%d" % (gai, gbi))
In [13]:
if DO_REGRESSIONS:
gbi, gai = (3, 0)
gB = groups_bio[gbi]
gA = groups_agg[gai]
nx, ny = len(gA), len(gB)
# ny, nx = len(gA), len(gB)
fig, axes = plt.subplots(ny, nx, figsize=(FSIZE*nx, FSIZE*ny),
sharex='col', sharey='row', gridspec_kw={'hspace': 0, 'wspace': 0})
regress_one_groups(gA, gB, A, head_agg, map_agg, B_log, head_bio, map_bio, collect, axes)
mk_plot("regression1VDent%dvClim%d" % (gai, gbi))
In [14]:
if DO_REGRESSIONS:
gbi, gai = (3, 1)
gB = groups_bio[gbi]
gA = groups_agg[gai]
nx, ny = len(gA), len(gB)
# ny, nx = len(gA), len(gB)
fig, axes = plt.subplots(ny, nx, figsize=(FSIZE*nx, FSIZE*ny),
sharex='col', sharey='row', gridspec_kw={'hspace': 0, 'wspace': 0})
regress_one_groups(gA, gB, A, head_agg, map_agg, B_log, head_bio, map_bio, collect, axes)
mk_plot("regression1VDent%dvClim%d" % (gai, gbi))
Multiple variables¶
In [15]:
if DO_REGRESSIONS:
fo = mk_out("regressionModels.txt")
for hb, xps in collect.items():
bi = map_bio[hb]
y = B_log[:, bi]
fo.write("\n# %s\n##############\n" % hb)
ks = sorted(xps.keys(), key=lambda x: regres_score(xps[x]))
fo.write("\n".join(["\t%s\t%s" % (v, regres_log(xps[v], [v], "")) for v in ks[:TOPR]])+"\n")
best_score = regres_score(xps[ks[0]])
new_xps = {}
new_compat = {}
for i, ki in enumerate(ks):
for j, kj in enumerate(ks[:i]):
pairH = (kj, ki)
pair = (map_agg[kj], map_agg[ki])
X = A[:, pair]
Xr = sm.add_constant(X)
# model = sm.OLS(y, Xr)
model = sm.GLS(y, Xr)
results = model.fit()
if regres_score(results) < best_score:
new_xps[pairH] = results
for vp in [0,1]:
if pairH[vp] not in new_compat:
new_compat[pairH[vp]] = set()
new_compat[pairH[vp]].add(pairH[1-vp])
while len(new_xps) > 0:
xps = new_xps
compat = new_compat
new_xps = {}
new_compat = {}
seen = set()
ks = sorted(xps.keys(), key=lambda x: regres_score(xps[x]))
fo.write("--------- (%d)\n" % len(ks))
fo.write("\n".join(["\t%s\t%s" % ("+".join(v), regres_log(xps[v], v, "")) for v in ks[:TOPR]])+"\n")
best_score = regres_score(xps[ks[0]])
for k in ks:
common = set.intersection(*[compat[kk] for kk in k])
for c in common:
vrsH = tuple(sorted([c]+list(k)))
if vrsH not in seen:
seen.add(vrsH)
vrs = [map_agg[v] for v in vrsH]
X = A[:, vrs]
Xr = sm.add_constant(X)
# model = sm.OLS(y, Xr)
model = sm.GLS(y, Xr)
results = model.fit()
if regres_score(results) < best_score:
new_xps[vrsH] = results
for vi, vp in enumerate(vrsH):
if vp not in new_compat:
new_compat[vp] = set()
new_compat[vp].update(vrsH[:vi])
new_compat[vp].update(vrsH[vi+1:])
close_out(fo)
Clustering¶
In [16]:
ctypes = ["C:ones", "C:sizes", "C:wdist"]
nkmeans = "k-means"
linkage_args = [("Ward", {"method": "ward"}),
("Complete", {"method": "complete"}),
# ("Single", {"method": "single"}),
# ("Average", {"method": "average"}),
# ("Centroid", {"method": "centroid"}),
("Weighted", {"method": "weighted"}),
("Median", {"method": "median"})]
collect_clusters = []
collect_clusters_names = []
linkZs = {}
Redescriptions based clusters¶
In [17]:
if DO_CLUSTERS:
for mi, ctype in enumerate(ctypes):
for ni, nbc in enumerate(NBCS):
keep_ids = C[:,map_clu["%s%d" % (ctype, nbc)]] > -1
cluster_labels = C[keep_ids, map_clu["%s%d" % (ctype, nbc)]].astype(int)
ccounts = numpy.bincount(cluster_labels)
collect_clusters.append(C[:, map_clu["%s%d" % (ctype, nbc)]].astype(int))
collect_clusters_names.append((ctype, nbc))
k-means clustering¶
In [18]:
if DO_CLUSTERS:
meths_args = [(nkmeans, cluster.KMeans, {"n_clusters": nbc, "random_state": 10}) for nbc in NBCS]
# ("AP", cluster.AffinityPropagation, {"max_iter": 50})]
for mi, (mname, meth, args) in enumerate(meths_args):
cluster_labels = meth(**args).fit_predict(clu_data)
ccounts = numpy.bincount(cluster_labels)[1:]
collect_clusters.append(cluster_labels)
collect_clusters_names.append((mname, args["n_clusters"]))
Hierarchical clustering¶
In [19]:
if DO_CLUSTERS:
for li, (lname, args) in enumerate(linkage_args):
Z = scipy.cluster.hierarchy.linkage(clu_data, **args)
linkZs[lname] = Z
for ni, nbc in enumerate(NBCS):
cluster_labels = scipy.cluster.hierarchy.fcluster(Z, nbc, criterion="maxclust")
ccounts = numpy.bincount(cluster_labels)[1:]
collect_clusters.append(cluster_labels)
collect_clusters_names.append((lname, nbc))
Ecoregions¶
In [20]:
if DO_CLUSTERS:
fig, axes = plt.subplots(figsize=(1.5*FSIZE, 1.5*FSIZE))
ecoregions_labels = C[:, map_clu["ECO_CODE"]].astype(int)
axes, bm, bm_args = tools_geomap.prepare_map(map_prefs, map_corners[0]+map_corners[1], axe=axes)
xs,ys = bm(lng_lat[:,0], lng_lat[:,1])
axes.scatter(xs, ys, c=ecoregions_labels, s=5, edgecolors='none', zorder=10, alpha=.7)#, cmap=cmC)
# axes.set_xlabel(" ".join(["%d" % len(ccounts)]+["%d" % i for i in ccounts]))
axes.set_title(f"Ecoregions")
mk_plot("clustersEcoregions")
Comparing clusterings¶
In [21]:
if DO_CLUSTERS:
nmi_meths = [ci for ci, nc in enumerate(collect_clusters_names) if nc[1] == NBC or nc[0] == nkmeans]
nmi_meth_names = ["%s %d" % collect_clusters_names[ci] for ci in nmi_meths]+["ecoregions"]
NMI = numpy.ones((len(nmi_meth_names), len(nmi_meth_names)))
for ii, ci in enumerate(nmi_meths):
for jj in range(ii):
NMI[ii,jj] = metrics.normalized_mutual_info_score(collect_clusters[nmi_meths[ii]], collect_clusters[nmi_meths[jj]])
NMI[jj,ii] = NMI[ii,jj]
NMI[ii,-1] = metrics.normalized_mutual_info_score(collect_clusters[nmi_meths[ii]], ecoregions_labels)
NMI[-1,ii] = NMI[ii,-1]
fig, axes = plt.subplots(1,1, figsize=(2*FSIZE, 2*FSIZE))
plt.imshow(NMI, cmap=cmB, vmin=-1, vmax=1)
plt.xticks(numpy.arange(len(nmi_meth_names)), nmi_meth_names, rotation=45)
plt.yticks(numpy.arange(len(nmi_meth_names)), nmi_meth_names)
axes.xaxis.set_ticks_position('top')
mk_plot("clustersNMI")
In [22]:
if DO_CLUSTERS:
clusters_sets = {}
map_same = {}
map_cc_cv = {}
clusterings_info = {}
for ci, cc in enumerate(collect_clusters):
ck = collect_clusters_names[ci]
clusterings_info[ck] = {"ci": ci,
"silhouette": metrics.silhouette_score(clu_data, cc)}
vals = set(cc)
map_cc_cv[ck] = sorted(vals)
for v in vals:
ss = set(numpy.where(cc==v)[0])
for okv, os in clusters_sets.items():
if okv[0] != ck and len(ss) == len(os) and ss == os:
if okv not in map_same:
map_same[okv] = okv
map_same[(ck, v)] = okv
break
if (ck, v) not in map_same:
clusters_sets[(ck, v)] = ss
cks = sorted(clusters_sets.keys(), key=lambda x: len(clusters_sets[x]))
contain = {}
for ci in range(len(cks)):
contain[cks[ci]] = []
for cj in range(ci):
if clusters_sets[cks[cj]].issubset(clusters_sets[cks[ci]]):
contain[cks[ci]].append((len(clusters_sets[cks[cj]])/len(clusters_sets[cks[ci]]), cks[cj]))
sub_parts = {}
for k, oss in contain.items():
if len(oss) > 1:
subs = {}
for fct, oks in oss:
if oks[0] not in subs:
subs[oks[0]] = []
subs[oks[0]].append(oks[1])
X = [okk for okk, ss in subs.items() if len(ss) > 1]
for x in X:
union_set = set().union(*[clusters_sets[(x, kk)] for kk in subs[x]])
if clusters_sets[k] == union_set:
if k not in sub_parts:
sub_parts[k] = []
sub_parts[k].append([(x, kk) for kk in subs[x]])
if k in sub_parts:
if len(sub_parts[k]) > 1:
print("Multiple mappings", k, sub_parts[k])
# pdb.set_trace()
sub_parts[k] = sub_parts[k][0]
sml_cks = [ck for ck in cks if ck not in sub_parts]
dists = numpy.zeros((len(sml_cks), len(sml_cks)))
for ci in range(len(sml_cks)):
for cj in range(ci):
dists[ci,cj] = 1-len(clusters_sets[sml_cks[ci]].intersection(clusters_sets[sml_cks[cj]]))/len(clusters_sets[sml_cks[ci]].union(clusters_sets[sml_cks[cj]]))
dists[cj,ci] = dists[ci,cj]
L = numpy.diag(dists.sum(axis=0)) - dists
egval, egvect = eigsh(L, 2)
ord_ks = numpy.argsort(egvect[:,1])#[-5:]
dist_seq = numpy.array([0.]+[dists[ord_ks[i-1],ord_ks[i]] for i in range(1, len(ord_ks))])
# dist_seq = numpy.array([0.]+[1*(dists[ord_ks[i-1],ord_ks[i]]>0) for i in range(1, len(ord_ks))])
scaleOD = numpy.cumsum(dist_seq/numpy.sum(dist_seq))
map_ccolor_scores = dict(zip([sml_cks[oi] for oi in ord_ks], scaleOD))
clusters_info = {}
for nc, cs in map_cc_cv.items():
for ci in cs:
ck = (nc, ci)
cck = map_same.get(ck, ck)
ccolor, parts = -1, []
if cck in map_ccolor_scores:
ccolor = map_ccolor_scores[cck]
elif cck in sub_parts:
parts = list(sub_parts[cck])
pp_in = [pk for pk in parts if map_same.get(pk, pk) in map_ccolor_scores]
pp_out = [pk for pk in parts if map_same.get(pk, pk) not in map_ccolor_scores]
while len(pp_out) > 0:
pparts = []
for pk in pp_out:
pparts.extend(sub_parts[map_same.get(pk, pk)])
pp_in.extend([pk for pk in pparts if map_same.get(pk, pk) in map_ccolor_scores])
pp_out = [pk for pk in pparts if map_same.get(pk, pk) not in map_ccolor_scores]
collect_ccolor = [(map_ccolor_scores[pk], len(clusters_sets[pk])) for pk in pp_in]
ccolor = sum([cc[0]*cc[1] for cc in collect_ccolor])/sum([cc[1] for cc in collect_ccolor])
parts = pp_in
else:
print("%s NOT FOUND!!" % ck)
clusters_info[ck] = {"basis": cck, "ccolor": ccolor, "parts": parts,
"center": numpy.mean(lng_lat[list(clusters_sets[cck]), :], axis=0),
"size": len(clusters_sets[cck]),
"label": chr(ord("A")+int(25*ccolor)) + chr(ord("a")+int(25*(25*ccolor)%25)),
}
In [23]:
if DO_CLUSTERS:
method_family = [nkmeans]+ctypes
nb_rows = len(NBCS)
nb_cols = len(method_family)
fig, axes = plt.subplots(nb_rows, nb_cols, figsize=(FSIZE*nb_cols, FSIZE*nb_rows))
for xi, nm in enumerate(method_family):
for yi, nbc in enumerate(NBCS):
clustering_info = {}
plot_clusters(axes, xi, yi, lng_lat, nm, nbc, map_cc_cv[(nm, nbc)], clusters_sets, clusters_info, clusterings_info[(nm, nbc)], cmC)
mk_plot("clustersKRMaps")
In [24]:
if DO_CLUSTERS:
method_family = [l[0] for l in linkage_args]
nb_rows = len(NBCS)+1
nb_cols = len(method_family)
fig, axes = plt.subplots(nb_rows, nb_cols, figsize=(FSIZE*nb_cols, FSIZE*nb_rows))
for xi, nm in enumerate(method_family):
scipy.cluster.hierarchy.dendrogram(linkZs[nm], p=3, orientation="left", truncate_mode="level", link_color_func=lambda k: "black", no_labels=True, ax=axes[-1,xi])
for yi, nbc in enumerate(NBCS):
clustering_info = {}
plot_clusters(axes, xi, yi, lng_lat, nm, nbc, map_cc_cv[(nm, nbc)], clusters_sets, clusters_info, clusterings_info[(nm, nbc)], cmC)
mk_plot("clustersHMaps")
